#ifndef CSPACECHARGE_H
#define CSPACECHARGE_H
#include <fstream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_halfcomplex.h>

#include <vector>
using namespace std;
const double clight=2.9929e8;
const double pi=3.14159;
class CSpaceCharge
{
    public:
	CSpaceCharge();
	~CSpaceCharge();
	//下面的函数可能不好实现。CTrackMp里面已经有beam，直接操作。这里就没办法把成员函数的beam当作整个0
	//程序的beam操作。如果用CTrackMp的beam初始化这个地方，可能会出现问题，初始化的是构造函数，但是
	//beam类的入口，不仅仅是构造函数，可能会出现问题。如果不使用beam的话，就要重新写关于rms尺寸、
	//
	void InputBeam(vector<vector<double> >*,double);
	void InGrid(int,int,int,double,double,double);
	void InputStruct(double _dFreq,double _dErest=938.2723);
	void setFreq(double _dFreq);
//	void setR(double _dRadius);
	void setCellLength(double _dLCell);
	void scheff(double);
   // void setFlag(int);
	void setPhis(double dPhis=0);
    void getVvv(vector<vector<vector<double> > > *vvv,string str,int iflag=0,double dTmpx=0,double dTmpy=0,double dTmpz=0,int m=0);
    //just for test:
   // void getBeam(int);
   // void getvvder();
   // void getvvdez();
   // void getComp(string);
   // void getPotenAlongR(string); //得到束团中心电场随着r的变化
    void formGrid(); 
    void FieldError();
   // void setStepNum(int);

    private:
	//void initField();
	void EFieldMap();
   // void flds(double, double, int, double, double, double, double*, double*);
   // void eint(double, double* , double*);
	void plim(double*, double* );
    void Interpret();
   // void EFieldMap2();
   // void Boundary(int);
	void Potential();
	void FFT3D(vector<vector<vector<double> > >*);
	void IFFT3D(vector<vector<vector<double> > >*);
	void FFT1D(vector<vector<vector<double> > >*);
	void IFFT1D(vector<vector<vector<double> > >*);
    vector<vector<vector<double> > >::iterator itervvvdBegin;
    vector<vector<vector<double> > >::iterator itervvvdEnd;
    vector<vector<double> >::iterator itervvdBegin;
    vector<vector<double> >::iterator itervvdEnd;
    vector<double>::iterator itervdBegin;
    vector<double>::iterator itervdEnd;
    public:
    //CBeam beam;
    //int stepNum;  //attention
    private:
    double xHalf,yHalf,zHalf;
	double dQ,dI;         //粒子电量和正负，现在考虑质子,dQ没用.dI is average current in one rf period.
//        double **compData;
    vector<vector<vector<double> > >compData,vvvdaa,vvvdex,vvvdey,vvvdez;
    int nx,ny,nz;
    double dx,dy,dz;
    vector<double> vdzm,vdxm,vdym,vdzs,vdys,vdxs;
	vector<vector<double> > vvdBeam;            //beam中束流粒子坐标分布
	vector<vector<double> > *cord;            //beam中束流粒子坐标分布
	double xbar,ybar,pbar;
	double dLCell,dPhis,dpLeft,dErest;        //cell同步相位和长度,同步能量
	double dFreq,dWaveLength;
	int iParticleNumTotal,ng;//nip表示计算相邻几个书团之间的空间电荷力，暂时仅考虑书团内的＝1
	double *dEnergyK,dBetaS,dGammaS;
    double dt;
    
    /*
	double dt,dz,zs,dr,hl,longigridlength,rmax,pl;            //网格尺寸，dt脉冲作用时间
	int nr,nz;              //网格信息
	int iflag;
	vector<double> vdrm,vdrs,vdrssq,vdzm,vdzzs,vders,vdezs; //网格坐标，场强,电荷加权值
	vector<vector<double> > vvdBeam;            //beam中束流粒子坐标分布
	vector<vector<double> > *cord;            //beam中束流粒子坐标分布
	vector<int>  vismax,viemax;
	double xbar,ybar,pbar;
	double *dEnergyK;
	double  dBetaS;
	vector<vector<double> > vvder,vvdez;                       //格点电量的加权值
	vector<vector<double> > vvdaa;
	int iBunchNum;                              //计算时考虑几个书团的作用
	double dLCell,dPhis,dpLeft,dErest,dRadius;        //cell同步相位和长度,同步能量
	double dFreq,dWaveLength;
        double ** dDensityK;
        double ** dDensity;
	int iParticleNumTotal,nip,np,ng;//nip表示计算相邻几个书团之间的空间电荷力，暂时仅考虑书团内的＝1
	//np 是计算过程中遍历用到的变量，ng表示好粒子束。不考虑书束流损失的时候等于iParticleNumTotal
	double c1,c2,c3,c4;
	double xOverR,yOverR,crp,czd;  //Interpret和scheff之间的数据传递
    */
};
#endif
